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We introduce a model for the spreading of epidemics by long-range infections and investigate the 
critical behaviour at the spreading transition. The model generalizes directed bond percolation and 
is characterized by a probability distribution for long-range infections which decays in d spatial 
dimensions as l/r d+a . Extensive numerical simulations are performed in order to determine the 
density exponent /3 and the correlation length exponents t/\i and v± for various values of a. We 
observe that these exponents vary continuously with a, in agreement with recent field-theoretic pre- 
dictions. We also study a model for pairwise annihilation of particles with algebraically distributed 
long-range interactions. 



PACS numbers: 05.70.Ln, 64.60.Ak, 64.60. Ht 
I. INTRODUCTION 

Spreading processes are often encountered in nature 
in situations as diverse as epidemics , catalytic reac- 
tions H, forest fires Q|, and transport in random me- 
dia HQ. Depending on the particular environmental 
conditions, the spreading process may either continue 
to spread over the whole population or die out after 
some time. The essential features of this transition be- 
tween survival and extinction of the spreading agent may 
be described by simple stochastic lattice models, which 
mimic the spreading mechanism by certain probabilistic 
rules. Usually such models incorporate two competing 
processes, namely spreading (infection) of nearest neigh- 
bors and spontaneous recovery (healing), with or with- 
out immunization. The spreading properties depend on 
the relative rates of the two processes. For example, if 
the rate for infection is very low, the spreading agent 
will disappear after some time and the system becomes 
trapped in an inactive state (or set of states) which is usu- 
ally referred to as the absorbing state of the model. On 
the other hand, if infections occur more frequently, the 
spreading process may survive for a very long time. The 
main theoretical interest in these models stems from the 
fact that the phase transition from the fluctuating active 
phase into the non-fluctuating absorbing state is contin- 
uous, and characterized by universal scaling laws associ- 
ated with certain critical exponents. As in equilibrium 
statistical mechanics, these exponents allow one to cat- 
egorize different lattice models into universality classes. 
Each of these universality classes then corresponds to a 
specific underlying field theory. 

The most important universality class for spreading 
transitions with short-range interactions is Directed Per- 
colation (DP) 0, as described by Reggeon field the- 
ory @ 0- in fact, DP covers the majority of phase 
transitions from a fluctuating active phase into absorb- 
ing states. The fundamental properties of DP have been 
expressed as a conjecture in Refs. pOlpll. Accordingly 



a spreading transition belongs to the DP class if (a) the 
absorbing state is unique, (b) the active phase is char- 
acterized by a one-component positive order parameter, 
(c) there are no other symmetries of the physical system 
except for spatio-temporal translation and spatial reflec- 
tion invariance, (d) there is no frozen disorder, and (e) the 
dynamical rules involve only short-range interactions. 

In many realistic spreading processes, however, short- 
range interactions do not appropriately describe the un- 
derlying transport mechanism. This situation emerges, 
for example, when an infectious disease is transported 
by insects. The motion of the insects is typically not a 
random walk, rather one observes occasional flights over 
long distances before the next infection occurs. Simi- 
lar phenomena are expected when the spreading agent 
is subjected to a turbulent flow. It is intuitively clear 
that occasional spreading over long distances will signif- 
icantly alter the spreading properties. On a theoretical 
level such a super-diffusive motion may be described by 
Levy flights i.e., by uncorrelated random moves over 
algebraically distributed distances. 

Anomalous directed percolation, as originally proposed 
by Mollison JjJ in the context of epidemic spreading, is 
a generalization of DP in which the spreading agent per- 
forms Levy flights. This means that the distribution of 
the spreading distance r is given by 

P(r) ~ l/r d+ ° , (a > 0) , (1) 

where d denotes the spatial dimension of the system. The 
exponent a is a free parameter that controls the charac- 
teristic shape of the distribution. It should be empha- 
sized that er does not introduce any new length scale, 
rather it changes the scaling properties of the underlying 
(anomalous) diffusion process. We are particularly inter- 
ested in the critical properties of anomalous DP close to 
the phase transition. As in the case of ordinary DP, we 
expect anomalous DP to be characterized by the univer- 
sal critical exponents (3, Vx, and u\\. The exponent 



1 



is related to the order parameter, the density of active 
sites n. Since the DP process is controlled by a single 
parameter p, with the phase transition taking place at 
p = p Cl then close to this transition in the active phase n 
vanishes as n ~ (p — p c f. At the same time, we expect 
the spatial and temporal correlation lengths £j_ and £|| 
to diverge as £_l ~ \p — p c \~ u± - and £|| ~ \p — p c \~ Uu , re- 
spectively, on both sides of the transition. Theoretically, 
one is interested in the dependence of these exponents 
on a, whether the exponents are independent from one 
another, and how they cross over to the exponents of or- 
dinary DP. Some time ago Grassberger jl2| claimed that 
the critical exponents of anomalous DP should depend 
continuously on the control exponent a. Very recently 
this work has been considerably clarified and extended 
by Janssen et. al. |r|| , who have presented a comprehen- 
sive field-theoretic renormalization group (RG) calcula- 
tion for anomalous spreading processes with and with- 
out immunization. The aim of the present work is to 
verify their results numerically. To this end we intro- 
duce a model for anomalous DP which generalizes di- 
rected bond percolation. In contrast to previously stud- 
ied models Jl4| , p^[ we do not introduce an upper cutoff 
for the flight distance r, and hence finite size effects are 
drastically reduced. Extensive numerical simulations are 
performed in order to determine the critical exponents, 
which are found to compare favourably with the field- 
theoretic predictions. 

The paper is organized as follows. In Section || wc 
first review the known field-theoretic results. In Sec- 
tion III we introduce a lattice model for anomalous DP 



and discuss the role of finite size effects. In Section |y| 
the numerical results are presented and compared with 
the field-theoretic predictions. We also discuss the case 
of anomalous pair annihilation in Section 



II. FIELD-THEORETIC PREDICTIONS 



In this section we will summarize some of the field- 
theoretic results which have been derived in Ref. [fl3|| . 
First of all let us recall that the Langevin equation for 
ordinary DP pOl is given by 



d_ 

at 



n(x,i) = (r + L>ArV 2 )n(x,t) - A n 2 (x, t) + C(x, t) , 

(2) 



where the constant r controls the balance between off- 
spring production and self-destruction, and plays the role 
of the deviation p — p c from the critical percolation prob- 
ability. The infection of nearest neighbors is represented 
by the diffusion operator V 2 , while the nonlinear term 
incorporates the exclusion principle on the lattice. The 



fluctuations are taken into account by adding a multi- 
plicative Gaussian noise field C( x i t) which is defined by 
the correlations 

(C(x,t)C(x',t')> =2r7i(x,i)<5 d (x-x')£(<-i') . (3) 

In order to generalize this Langevin equation to the case 
of anomalous DP, the short-range diffusion has to be re- 
placed by a non-local integral expression which describes 
long-range spreading according to the probability distri- 
bution P(r): 



d_ 

dt 



t(x, t) = t n(x, t) — A n 2 (x, t) + £(x, i) 



(4) 



+D /^P( |x -x1)Kx', t )-„(x, t )]. 



The two contributions in the integrand describe gain and 
loss processes, respectively. Keeping the most relevant 
terms in a small momentum expansion, this equation may 
be written as [|13| 



-|n(x, t) = (D N V 2 + D A V° + r ) v ( x. / i 
-An 2 (x,t) +C(x,i) , 



(5) 



where the noise correlations are assumed to be the same 
as in Eq. ([|). Djy and Da arc the rates for normal and 
anomalous diffusion, respectively. The anomalous diffu- 
sion operator V CT describes moves over long distances and 
is defined through its action in momentum space 



V 7 . 



(6) 



where k = |k|. The standard diffusive term -DatV 2 takes 
into account the short range component of the Levy dis- 
tribution. Note that even if this term were not initially 
included, it would still be generated under renormaliza- 
tion of the theory. 

Before summarizing the field-theoretic results, let us 
first consider the mean-field approximation. As in ordi- 
nary DP the mean-field dynamic phase transition occurs 
at r = 0, where gain and loss processes balance one an- 
other. For t < 0, the particle density decays exponen- 
tially quickly towards n — 0, which is the absorbing state 
of the system. However, for t > 0, the stable stationary 
state now has the non-zero particle density n = r/A. 
Since r plays the role of p — p c , the mean field density 
exponent is fj MF = 1. The scaling exponents u± and f|| 
can be derived from an inspection of Eq. (H). For a < 2, 
we see that 
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,MF 



= 1/<T , 



(7) 



and the characteristic time diverges according to 



2 



,MF 



(8) 



As expected, for a > 2, these exponents cross over 
smoothly to the ordinary DP exponents. Note that the 
mean field result demonstrates that u± varies continu- 
ously with a . 

The mean field approximation is expected to be quan- 
titatively accurate above the upper critical dimension. 
For d < d c , however, fluctuation effects have to be taken 
into account. The fluctuation corrections to the critical 
exponents can be computed by a field-theoretic RG cal- 
culation. Using standard techniques, the Langevin equa- 
tion (|^) can be rewritten as an effective action: 



dt 



$(d t - t - D N V 2 - D A V a )^ 
+ |(^ 2 -^)]- (9) 



Simple power counting on this action reveals that the 
upper critical dimension is d c — 2cr, below which fluc- 
tuation effects become important. In the above action 
the field ^(x, t) can be identified with the coarse-grained 
particle density field n(x, t) |]l6| and ^(x, i) is the corre- 
sponding response field. The expression in Eq. (Mj differs 
from the usual action of Reggeon field theory p|-[To|] by 
the addition of a term representing anomalous diffusion. 

The field-theoretic RG calculation in Ref. jl3| em- 
ploys Wilson's momentum shell renormalization group 
recursion relations in order to determine the critical ex- 
ponents. The authors of the present work have inde- 
pendently performed similar calculations based on di- 
mensional regularization which are fully consistent with 
Ref. [fl3[ . In the following we summarize the main results. 
The critical exponents to one-loop order in d = 2a — e 
dimensions are given by 

/3 = l-| + 0(> 2 ) , 

^ = I + il + o( e 2 ) , 

a la £ 



i + ^ + o( e 2 ) 

a-e/7 + 0(e T 



(10) 



Moreover, it can be shown that the hyperscaling relation 
9 + 2S = d/z (<5 = /3/V||) , (11) 

for the so-called critical initial slip exponent 9 |l7j , holds 
for arbitrary values of a. The exponent 9 (which is also 
sometimes denoted by r] in the literature) describes the 
initial increase in the number of active particles N(t) for 
critical systems starting from initial states at very low 



density, i.e. where we have N(t) ~ t B . The critical initial 
slip plays an important role in dynamical Monte-Carlo 
simulations (see Section [Tv|) . To one- loop order, 9 and 5 
are given by 



-£+<>(«■) 



«=1-§ + 0(eV (12) 



Finally, thanks to the fact that Da does not get renor- 
malized, one can prove the exact scaling relation 



/|| - uj_(a - d) - 2/3 = 



(13) 



The fact that Da does not get renormalized means that 
anomalous DP is described by two rather than three in- 
dependent critical exponents. The scaling relation ( |l3| ) 
has a further surprising consequence. Assuming that /?, 
i>j_ and v\\ change continuously with <r, then for fixed d, it 
predicts the value a c where the system should cross over 
to ordinary DP (assuming the crossover is smooth). To 
this end one simply has to insert the numerically known 
values of the DP exponents into Eq. ([l3|). Surprisingly 
one obtains a c = 2.0766(2) in one, a c ~ 2.2 in two, and 
cr c = 2 + e/12 in d = 4 — e spatial dimensions. Thus the 
crossover takes place at a c > 2 which collides with the 
intuitive argument that the anomalous diffusion opera- 
tor V 17 should only be relevant if a < 2. But, as pointed 
out in Ref. |]l3| , this naive argument may be wrong in 
an interacting theory where the critical behaviour is de- 
termined by a nontrivial fixed point of an RG transfor- 
mation. Rather the field-theoretic calculation predicts 
that anomalous diffusion is still relevant in the range 
2 < a < <r c (d) for d < 4. This prediction seems to 
be additionally surprising since the operators for anoma- 
lous and ordinary diffusion V CT and V 2 are expected to 
coincide for a — 2. However, one can show that for a = 2 
the most relevant terms in a small momentum expansion 
of Eq. M) also contain a logarithmic correction of the 
form — Arlogfc. Therefore anomalous and ordinary dif- 
fusion are indeed different in that case, supporting the 
view that long-range spreading might be relevant in the 
regime 2 < a < a c . Unfortunately, the numerical simu- 
lations presented in Section ^ are not accurate enough 
to confirm this prediction. 

Another interesting aspect of anomalous DP is that 
o~ can be chosen in such a way that the critical dimen- 
sion d c = 2a approaches the actual physical dimension 
at which the simulations are performed (see Section |^) . 
Even in one spatial dimension this allows us to verify 
the one-loop results (|l(]). For example, if a = 1/2 + fj,, 
the critical dimension of the system is d c = 1 + 2[i and 
hence the exponents in a 1+1-dimcnsional system change 
to first order in /i as 

(3 = 1 - 8^/7 + O {y 2 ) , 
v ± = 2 - 12/i/7 + O ( M 2 ) , 
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(a) Ordinary directed bond percolation: 



t 

t+1 



• • • • A* ' " * 
.......... 



(b) Anomalous directed bond percolation: 
t • • • • • • • • • • • 

t+1 . y ....... \ 



2d L -l 



2d R -1 



FIG. 1. Dynamical rules for (a) ordinary directed bond 
percolation and (b) the present model with algebraically dis- 
tributed distances. d,L and cLr are defined in the text. 



v\\ = 1 + 4/i/7 + O (m 2 ) , 
z = 1/2 + 5a*/7 + O (a* 2 ) 
5 = 1- 12^/7 +0(fx 2 ) 
9 = 4^/7 + O (fi 2 ) . 



(14) 



In Section IV we shall demonstrate that this initial 



change of the exponents is indeed observed in numeri- 
cal simulations. 



III. A LATTICE MODEL FOR ANOMALOUS 
DIRECTED PERCOLATION 



Anomalous DP was first studied numerically by Al- 
bano |l4|| who introduced a model for branching-annihi- 
lating random walks in which the particles performed 
Levy flights. However, his estimates for the critical ex- 
ponents were rather inconclusive, in particular they vi- 
olated the scaling relation (|l3| ) and even the mean-field 
limit was not correctly reproduced. In Ref. (l3| it was 
suspected that these problems could have originated in 
the truncation of the flight distances at some upper cut- 
off, usually at the system size. The upper cutoff effec- 
tively suppressed long range motion and hence DP-likc 
behaviour was amplified. A systematic finite size analysis 
of a similar model confirms this point of view and shows 
that even on a lattice with 10 4 sites finite-size effects are 
still extremely dominant. 

Similar problems were also encountered in a more re- 
cent study of a generalized Domany-Kinzel model with 
long range interactions (lq |. In this case an upper cutoff 
for P(r) was also introduced (by defining transition prob- 
abilities w(5*|S'*~ 1 ) in which the sum over the spreading 
distance for a system with N sites is truncated at N/2). 
It was reported that the percolation threshold depended 
on the system size and varied by more than 20%. How- 
ever, it seems that this unusual drift of p c is actually 





a=°o 



a=2.0 



v.^:. ^' .* * •.• - , 



o=1.0 a=0.5 

FIG. 2. Critical anomalous directed percolation in 1+1 
dimensions for different values of a. The figure shows typical 
clusters starting from 5 active sites in the center of the lattice. 
The case a — oo corresponds to ordinary DP. As a decreases, 
spatial structures become more and more smeared out until 
in the mean field limit a = 1/2 the particles appear to be 
randomly distributed over the whole system. For small values 
of a finite size effects may lead to sudden transitions into the 
absorbing state. 

related to extremely strong finite size effects. 

In order to minimize finite size effects, we introduce 
a model in which the probability distribution for long- 
range spreading is not truncated from above. As in the 
case of ordinary directed bond percolation, our model is 
defined on a tilted square lattice and evolves by paral- 
lel updates. A binary variable Si(t) is attached to each 
lattice site i. Sj = 1 means that the site is active (in- 
fected) whereas Sj = denotes an inactive (healthy) site. 
Although the model may be defined in arbitrary spatial 
dimensions, we will focus here on the 1+1-dimensional 
case. The dynamical rules (see Fig. [j]) depend on two 
parameters, namely the control exponent a > and the 
bond probability < p < 1. For a given configuration 
{si(t)} at time t, the next configuration {si(t + 1)} is 
constructed as follows. First the new configuration is 
initialized by setting Si(t + 1) := 0. Then a loop over all 
active sites i in the previous configuration is executed. 
This loop consists of the following steps: 

1. Generate two random numbers zl and zr from a 
flat distribution between and 1. 
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2. Define two real-valued spreading distances tl = 
z L 1 ' a and tr = z^'° , for spreading to the left (L) 
and to the right (R). The corresponding integer 
spreading distances cLl and are defined as the 
largest integer numbers that are smaller than 
and rn, respectively. If cLl or da exceed the allowed 
range for integer numbers we go back to step 1. 

3. Generate two further random numbers ul and yn 
drawn from a flat distribution between and 1, 
and assign s i+ i- 2 d L (t + 1) := 1 if tjl < P, and 
Si-i+2d R (t + 1) := 1 if Vr < p, respectively. In 
finite systems the arithmetic operations in the in- 
dices are carried out modulo L by assuming peri- 
odic boundary conditions, i.e. s< = Si±L- 

The model includes two special cases. For a — > oo 
it reduces to ordinary directed bond percolation with 
p c ~ 0.6447. On the other hand, for a — » the in- 
teraction becomes totally random. In that case the 
model is exactly solvable and the transition takes place 
at p c = 1/2. In between, the spreading properties of the 
model change drastically, as illustrated in Fig. 0. 

As can be easily verified, the assignment r — z^ 1 ! 
reproduces the normalized probability distribution 

f a/r 1+a if r > 1 , 
| otherwise . 

As usual the distribution has a lower cutoff at r m i n = 1, 
which represents the lattice spacing. But in contrast to 
previously studied models, no upper cutoff is introduced 
and therefore almost arbitrarily large spreading distances 
may be generated (limited only by the maximal range 
of 64-bit integer numbers). In finite systems the tar- 
get site is determined by assuming periodic boundary 
conditions, i.e., the particle may "revolve" several times 
around the system. It turns out that this simple trick 
considerably reduces finite size effects. In particular the 
value of p c is well defined over a wide range of system 
sizes. Nevertheless finite size effects are still important 
in this model. In particular for small values of a, where 
long-distance flights occur frequently, finite-size effects 
enhance the probability for a target site to be already 
occupied. This in turn reduces the average density of ac- 
tive sites in a growing cluster and therefore increases the 
probability to enter the absorbing state. For example, 
for a = 0.5, a small system with only 200 sites reaches 
the absorbing state typically after only a few hundred 
time steps (see Fig. ||). Therefore, in order to further 
reduce finite size effects, we either use very large lattice 
sizes of about 10 5 sites (as in stationary simulations, see 
below) or else in other (dynamical) simulations, we can 
eliminate finite size effects almost completely by working 
on a virtually infinite lattice (see also below). 




FIG. 3. The critical percolation threshold of anomalous 
directed bond percolation as a function of a. 

IV. NUMERICAL RESULTS 

In order to estimate the critical exponents of anoma- 
lous DP we employ two different standard Monte-Carlo 
techniques, namely dynamical simulations at criticality 
and steady-state simulations in the active phase. 

In dynamical simulations ]l8| a critical cluster is grown 
from a single active seed (just as in Fig. Averag- 
ing over many independent realizations one measures the 
survival probability P(t), the number of active particles 
N(t), and the mean square spreading of surviving clus- 
ters from the origin R 2 (t) . At criticality, these quantities 
are expected to scale as 

P(t)~t~ 5 , N{t)~t 9 , R 2 {t)~t 2 / Z , (16) 

where 5 = @fv\\ and 9 is the critical initial slip expo- 
nent |M (see Eq. ([Tl|)). Since the size of the growing 
cluster is finite, we are able to perform the simulations 
on a virtually infinite lattice by storing the coordinates 
of active particles in a dynamically generated list. The 
effective system size is then determined by the maximal 
spreading range (i.e., the maximal range of integer num- 
bers ±2 63 ), which means that finite size effects are al- 
most eliminated. Since deviations from criticality lead 
to a curvature of P(t) in a double logarithmic plot, the 
dynamical simulation method allows a precise estimate 
of the percolation threshold p c for different values of a 
(see Fig. || and Table |). As expected, p c tends to 1/2 in 
the limit a — > 0. 

Having determined the critical points, we measure the 
quantities P(t), N(t), and R 2 (t) at criticality. However, 
it turns out that, in the presence of sufficiently long-range 
interactions, the mean square spreading, defined as an 
arithmetic average R 2 (t) = (|x(i)| 2 ), diverges. In order 
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0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 

o a 



FIG. 4. Estimates for the critical exponents from dy- 
namical Monte-Carlo simulations in comparison with the 
field-theoretic predictions (solid lines) and the DP exponents 
(dot-dashed lines). 

to circumvent this problem, we instead compute the ge- 
ometric average 

i? 2 (t) = exp[(log(|xW| 2 ))] . (17) 

This average turns out to be finite for all a > and 
renders consistent results in the case of ordinary DP. The 
numerical estimates for the dynamical exponents 8, 8, 
and z/2 are shown in Fig. || 

The exponent (3 is determined by stationary simula- 
tions in the active phase. As the active phase of anoma- 
lous DP is characterized by a homogeneous particle den- 
sity, this type of simulation has to be performed on a 
finite lattice. In order to minimize finite size effects, we 
choose a large lattice size of L — 10 5 sites. Starting from 
a fully occupied initial state, the system first equilibrates 
over 10 time steps before the stationary density n is av- 
eraged over another 10 4 time steps. Our estimates for j3 
are shown in Fig. |^. Combining the results we can now 
compute the scaling exponents v± = /3/Sz and v\\ — (3/5, 
which are summarized in Table |. 

According to Eq. (p^), the one-loop expansion pre- 
dicts the initial variation of the critical exponents close 
to a = 1/2. This is one of the rare cases where one 
can directly "see" the field-theoretic results in the simu- 
lation data. In Figs. [H| the predicted initial slopes are 
indicated by solid lines. Clearly they are in fair agree- 
ment with the numerical estimates, which confirms the 
field-theoretic results of Ref. pa]. For a > 1.5, however, 
the numerical results are not accurate enough to verify 
the predicted location of the crossover to ordinary DP 
at a c — 2.0766(2). It seems that the deviations in this 



FIG. 5. Estimates for the exponent f3 and the derived 
exponents u± and v\\ in comparison with the field-theoretic 
results (solid lines) and the DP exponents (dot-dashed lines). 
The quantities Ai and A2 represent deviations from the scal- 
ing relations (|ll| and (^), respectively (see text). 

regime are due to very long crossover times in the dynam- 
ical simulations, resulting from a complicated interplay 
between long-range and short-range processes. 

In order to verify the scaling relations (|il| ) and ( fl3| ) 
we have also plotted the deviations Ai = 2£ + # — 1/z 
and A2 = 1 — <t+ (1 — 25) z which should be equal to zero 
in the intervals a > 0.5 and 0.5 < a < o~ c , respectively. 
In fact, as shown in Fig. [5], Ai is smaller than the error 
tolerance, which confirms the validity of the hyperscal- 
ing relation (|ll|). Similarly the values of A2 confirm the 
validity of Eq. ( |l3| ) in the range 0.5 < a < 1.5, whereas 
significant deviations occur for a > 1.5. We believe that 
these deviations do not indicate that the scaling rela- 
tion (|l3]) is violated for large values of <r, rather they 
confirm that the simulations in this regime may be af- 
fected by very long crossover times. 



V. ANOMALOUS ANNIHILATION PROCESS 

In this section we consider the somewhat simpler case 
of anomalous pair annihilation A + A — > with long- 
range hopping. This model was previously studied in 
p9[ , using both simulations and approximate theoretical 
techniques. In this paper we will extend this previous 
work, by presenting a systematic field-theoretic analysis, 
as well as by performing more detailed numerical simu- 
lations. 

In the ordinary annihilation process |20| with short- 



6 



a 


Pc 





v±_ 




2 


e 


0.2 


0.50026(1) 


0.99(4) 


4.7(6) 


1.00(6) 


0.21(2) 


0.01(3) 


0.3 


0.50097(1) 


0.98(5) 


3.2(3) 


1.01(6) 


0.31(2) 


0.02(3) 


0.4 


0.50245(2) 


0.97(6) 


2.5(2) 


0.99(7) 


0.39(2) 


0.01(3) 


0.5 


0.50490(2) 


0.95(6) 


1.87(13) 


1.01(6) 


0.54(2) 


0.02(3) 


0.6 


0.50847(2) 


0.88(5) 


1.76(12) 


1.02(6) 


0.58(2) 


0.04(3) 


0.8 


0.51820(3) 


0.76(4) 


1.60(11) 


1.13(7) 


0.71(2) 


0.13(3) 


1.0 


0.52981(5) 


0.65(3) 


1.52(11) 


1.25(7) 


0.82(2) 


0.20(3) 


1.2 


0.54197(5) 


0.56(3) 


1.46(10) 


1.40(9) 


0.96(3) 


0.28(2) 


1.4 


0.55390(10) 


0.49(3) 


1.36(11) 


1.48(11) 


1.09(3) 


0.30(2) 


1.6 


0.56520(10) 


0.43(3) 


1.29(11) 


1.56(14) 


1.21(3) 


0.30(2) 


1.8 


57561 < W) 


0.39(3) 


1 23 ( IS) 


1 62f 18^1 


1.32(3) 


0.31(2) 


2.0 


0.58505(10) 


0.34(3) 


1.13(15) 


1.62(20) 


1.43(3) 


0.32(2) 


2.2 


0.59345(10) 


0.32(3) 


1.13(15) 


1.68(21) 


1.49(4) 


0.31(2) 


2.4 


0.60085(10) 


0.30(3) 


1.15(17) 


1.76(24) 


1.53(5) 


0.32(2) 


DP 


0.644700 


0.2765 


1.097 


1.734 


1.581 


0.3137 



TABLE I. Estimates of the percolation threshold and the 
critical exponents for various values of a, compared to the 
corresponding values for ordinary bond DP. 



range interactions, the average particle density is known 
to decay as 



n(i) 




t -d/2 for d < 2 5 

"Mni for d = d c = 2 
- 1 for d > 2 . 



(18) 



Hence, except for the log correction in d = 2, the den- 
sity decays away as a power law, n(t) ~ t~ a . Turn- 
ing now to the Levy-flight case, this may be described 
theoretically by inserting an additional operator V 7 into 
the well-known field-theoretic action for pair annihilation 
(see pfj[). The resulting action reads 



S$,*l>] 



d d xdt 



{i>{d t - D N V 2 - DaW)^ 
+2X^j 2 + A^V - n #(i)} , (19) 



where uq is the initial (homogeneous) density at t = 0. 
Here the field ip is not simply related to the coarse- 
grained density field although it is true that the 
average values of both fields are the same. The action 
( |l9| ) can be derived systematically, starting with an ap- 
propriate (non-local) Master equation — the details are 
given in Appendix A. Note also that the action for the 
process A + A — > A with Levy flight hops differs only 
in the coefficients of the reaction terms. Hence the Levy 
flight annihilation and coagulation processes are in the 
same universality class. 

An analysis of the above action follows very closely 
that of Ref. |20| ], For a < 2, power counting reveals 
that the upper critical dimension of the model is now 
d c = a < 2. For d > d c mean- field theory is expected 
to be quantitatively accurate, with an asymptotic den- 
sity decay ~ t . Below the upper critical dimension, 



however, the renormalized reaction rate flows to an or- 
der e = <7 — d fixed point. This allows us to very quickly 
determine the asymptotic density decay via dimensional 
arguments. Below d c the only dimensionful quantity left 
in the problem is the time t, which, for a < 2, scales as 
[t] ~ fc _<T . Hence, for a < 2, the density must decay as: 



t- d / a for d < a , 
n(t) ~ { t~ l hit for d=d c 
t" 1 for d > a . 



(20) 



The derivation of the logarithm at the upper critical di- 
mension requires a slightly more sophisticated calcula- 
tion, which is, however, completely analogous to that in 
Ref. f2(i[| . Note also that for a > 2 the results cross 
over smoothly to the standard annihilation exponents of 
Eq. ©. 

A lattice model for anomalous annihilation in 1 + 1 
dimensions may be constructed by a simple modification 



of the model for anomalous DP introduced in Section III . 
To this end steps 1-3 have been modified such that for 
all active sites i we perform the following procedure: 

1. Generate a random number z 6 [0, 1] and define a 
real- valued spreading distance r = z~ x l a ' . The cor- 
responding integer spreading distance d is defined 
as the largest integer number smaller than r. 

2. Generate another random number y £ [0, 1] and 



assign «^_ M (t+l) := 



and 



-l+2<2 



1) := I - sfA^ 4 {t 



„old 
d 

-l+2d 



M (t+l)if»>l/2, 
1) otherwise. 



As in the case of anomalous DP, the arithmetic 
operations in the indices are carried out modulo 
L by assuming periodic boundary conditions, i.e. 
Si = s l±L . 

In step 2 the state of the target site is simply inverted. 
Therefore, if two particles move to the same target site, 
they annihilate instantaneously. 

Since the annihilation process starts with a homoge- 
neous density of particles, it is impracticable to work 
with a virtually infinite lattice by storing the coordinates 
of individual particles in a list. Rather we have to per- 
form ordinary simulations with a fixed system size. In 
order to minimize finite size effects we choose a large lat- 
tice size of 2 16 sites. For various values of a we measure 
the particle density n(t) up to 10 4 time steps averaged 
over at least 10 3 independent runs. By measuring the 
slopes of n(t) in a double logarithmic plot in the decade 
10 3 < t < 10 4 , we estimate the density decay exponents 
a(cr), which are shown in Fig. ^ (labeled as direct mea- 
surements). For a < 1.5 the agreement with the theoret- 
ical result of Eq. ( pC| ) (the solid line) is quite convincing, 
whereas large deviations occur close to a = 2. A de- 
tailed analysis of the local slope of n(t) as a function of 
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0.4 1 1 1 1 1 1 1 1 ' 1 1 

0.5 1 1.5 2 2.5 3 

a 



FIG. 6. The anomalous annihilation process: the graph 
shows direct estimates and extrapolations for the decay expo- 
nent a, as a function of a. The solid line represents the exact 
result (neglecting log corrections at a = 1). 

time in a double-logarithmic plot shows that these devi- 
ations are related to very long crossover times. In fact, 
determining the local slopes of n(t) in a log-log plot and 
extrapolating them graphically to t — > oo, one obtains a 
much better coincidence. On the other hand, for a < 1, 
the extrapolation leads to larger deviations. These er- 
rors may be related to finite size effects which are still 
extremely dominant in this regime. 



VI. CONCLUSIONS 



law waiting time distribution for the particles, in addi- 
tion to the power law Levy distribution for particle hops. 
Hence, in the absence of interactions, each particle would 
perform a continuous time Levy flight (see pl| and refer- 
ences therein). This modification should lead to a further 
universality class, with the exponents depending contin- 
uously on the control parameters for both the Levy-flight 
and the (power law) waiting time distributions. 



APPENDIX A: DERIVATION OF THE 
ANOMALOUS ANNIHILATION ACTION 

In this appendix we briefly describe how the anoma- 
lous annihilation action ( fl9| ) may be derived. To simplify 
matters we will not include the reaction terms — their 
derivation is precisely the same as in p0[ . The appro- 
priate Master equation for (anomalous) diffusion is given 
by 

|p(M; *) = ^£E [fa + i) qji (Ai) 

xP(. . .,rij + 1,. .. ,rii - 1, .. .;t) - riiqij P({n};t)] , 

where P({n}) is the probability of the particle configu- 
ration {n} = (rii, . . . , njv), I is the microscopic lattice 
spacing, Do is the diffusion constant, and where qtj gives 
the appropriate weight for a hop from site i to site j. Fol- 
lowing p^| , we next introduce creation and annihilation 
operators aj, a\, such that 

a\\rii) = \rii + 1) , ai\i%i) = m\rii - 1) , (A2) 



In this paper we have, for the first time, studied numer- 
ically the behaviour of Levy-flight DP close to the phase 
transition between the active and absorbing states. To 
this end, we have introduced a lattice model for anoma- 
lous DP in which finite size effects are considerably re- 
duced. Using special simulation techniques we have ob- 
tained accurate values for the associated critical expo- 
nents in 1 + 1 dimensions, which are almost free from 
finite-size effects. In addition we have performed quanti- 
tative tests of the one loop field-theoretic results, by tun- 
ing the upper critical dimension of the model to lie just 
above d = 1 . Our results are all in good agreement with 
the recent field-theoretic analysis of jlj|. Close to a = 2, 
however, our numerical results are not accurate enough to 
confirm the form of the predicted crossover to ordinary 
DP. We have also considered the simpler case of Levy- 
flight pair annihilation, where our numerics are again in 
agreement with (exact) field-theoretic arguments. 

Various possible extensions of the above models are 
possible. The most obvious involves including a power 



with the commutator [a,-,at] = Sij. The system state is 
then given by 

!*(*)>= E nw;*)*r •••<*£>>. (as) 

ni ,. . . ,njv 

where |0) is the vacuum state. Hence we can rewrite 
Eq. @) as 



where 



d_ 

di 



]a i. i 



*(*)> = -n\m) 



(A4) 



(A5) 



We may now perform the mapping to a field theory using 
standard methods (see |22| ) . After taking the continuum 
limit in space, we end up with the continuum action 
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S = J d d x dt {^(x, t) d t ip(x, t) - Di ^{x, t) J d d y 

x{[^{y,t)-^(x,t)}f\x-y\)} , (A6) 

where we have the Levy distribution 



f(r) d d r 



1 



r (T-\-d 



d d 7 



Transforming this into Fourier space, we obtain 
d d k 



(2tt 



- t dt { V(fc, t)dti){-k, t) - D x [iP(k, t) 

x(f(k)-l)i>(-k,t)}} , (A8) 



with 



f(k) - 1 = — I d d r 



r d-\-cr 



(A9) 



where AT is a normalization constant. After some ma- 
nipulation of the above integral, and after performing a 
small momentum expansion, we end up with 



S = 



d d k 
W) 



2 dt [4>(k,t)dS(-k,t) 



+4>(k, t){D A k a + D N k 2 + 0(A: 4 )}^(-fc, *)] , (A10) 

valid for < a < 4, a ^ 2. The final action ( p~9j ) is then 
obtained by the inclusion of both the reaction terms, and 
the initial density source. 
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